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Abstract 

In this paper we analyse the drying of a soil composed of particles, 
water and solute impurities, and study the occurrence of convective in- 
stabilities during evaporation. We find that the main driving force for 
instability is the formation of a concentration gradient at the soil surface 
due to the evaporation of water. A similar phenomenon may occur during 
the thawing of frozen ground in Arctic regions. 

1 Introduction 

The formation of a surface seal in soil, the hardening of the seal to form a crust 
and cracking of this crust, has been observed to occur under the combined in- 
fluence of rain and dry weather [TJ [2] . The sealing leads to a local increase in 
the bulk density and a decrease in porosity along with a decrease in the hy- 
draulic conductivity [3| . After drying these surface seals transform into surface 
crusts. The sealing and crusting of soil surfaces increases the runoff in the soil 
surfaces. In fact, initially the soil loss reduces because soil strength increases 
during crusting. However, crusting eventually leads to the formation of cracks 
which increases runoff in the soil surfaces [4]. The detailed effect of crusting 
on soil loss depends on various factors like the distribution of crusted and non- 
crusted regions and the shrinkage cracks within the crusted layer [5], Cracking 
depends critically on the rate of drying and so it has been observed to occur 
within a few days of a rainstorm in a dry season [5] . 

The effects of sealing and crusting on soil loss are most visible in agricultural 
lands. This is because the sealing and crusting is most likely to occur in the soil 
which is not vegetated in the agricultural environment. In Europe a quarter of 
its agricultural land exhibits some form of soil erosion risk [5]. In fact, a 20 mm 
rain has been observed to form surface seal and cause soil loss in certain hilly 
areas [9J [10] . In this case, the thickness of the surface seal was observed to be on 
the order of a few millimeters. However, the formation of a 7 cm thick surface 
seal during a wet season has also been observed [TT] . Since the soil loss can effect 
the properties of soil and thus have a major impact on agricultural production, it 
is important to understand the physical processes at work in drying ground. In 
this paper we analyse convective instabilities that are observed to occur during 
the evaporation of a mixture of soil particles and water [22] . 
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Evaporation of pure liquids can cause convection in those liquids [T3I [T3J 
RBI IT51 [TCI Ul] . This is because the evaporation causes the surface tempera- 
ture to drop and thus sets up an unstable density gradient which can induce 
a Rayleigh-Bernard instability. In addition, for most liquids surface tension is 
a function of temperature and so perturbations in temperature along the film 
surface also create perturbations in the surface tension. The liquid then flows 
from places of lower surface tension to places of higher surface tension. This 
flow is called the Marangoni effect. Thus, both the temperature gradient and 
surface tension variations contribute to the occurrence of convection in pure 
liquids. The combined phenomenon is called Bcrnard-Marangoni convection. 

The formation of surface seals and crusts can be understood by analysing the 
drying of pools of muddy water. Evaporation in mixtures also leads to convec- 
tion JTH1 HH1 HOI [3J • In this case as the liquid drys up a concentration gradient 
is also set up along with a temperature gradient. The liquid moves under the 
combined influence of both these gradients. Furthermore, the Marangoni effect 
also occurs due to the dependence of the of surface tension on both the temper- 
ature and the concentration. Thus, a perturbation in either the temperature or 
the concentration can cause convection. Convection has been observed to occur 
in the course of evaporation of a binary mixture of soil particles and water [33] . 
In doing so a four weight percent mixture of bentonite and water was analysed 
at twenty degrees centigrade. It was found that convection occurred during the 
course of drying of this mixture and produced hexagonal patterns. However, it 
was also observed that the temperature gradient set up during the evaporation 
could not explain the occurrence of convection in this system. In this paper we 
will show that the concentration gradient is the main driving force behind such 
convection during the drying of ground. 

We also study the thawing of frozen ground which occurs in the Arctic circle. 
Freezing and thawing of soil in the Arctic circle results in the formation of 
various surface patterns such as soil hummocks and stone circles [33] 12H HiJ [315] . 
Various models have been proposed to explain the occurrence of these patterns. 
In one of these models the maximum density of water at four degree centigrade 
sets up the convection [27] [28] [29]. This model is based on the convection 
of water through soil pores which has not been observed for the soils under 
consideration. It has also been proposed that the sedimentation of soil during 
thawing sets up an unstable density profile which causes convection [3D] [3T] [33] . 
However, the measured soil density gradient is not large enough to initiate 
convection |33j . Another model proposes that the motion responsible for pattern 
formation occurs during the freezing of ground (differential frost heave) [34j [35] . 
This model yields a plausible mechanism for patterned ground but has yet to 
be experimentally confirmed, and does not explain certain observations such as 
the soil convection observed to occur in later summer |25j . 

The soil convection can potentially be explained by a phenomenon similar 
to the one that occurs in the formation of surface seals. The ice melts in late 
summer and at the same time evaporation takes place from the surface of the 
soil. This causes an unstable density profile to develop that is in principle 
large enough to initiate wholesale soil convection. This convection has been 
observed in fields and is thought to contribute to the patterns that form in 
Arctic region [25] [33] . It may be noted that similar patterns have been detected 
on Mars, suggesting that in the past the temperature on Mars may have been 
large enough for the ice to melt and soil convection to take place [36) . 
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In this paper we will use tensor notations for performing the stability analysis 
[57] . We will also use a new approach based on Green's functions to perform 
stability analysis. As will be evident that the stability analysis of a ternary 
mixture is very complicated, and it would not be possible to perform it using 
the conventional methods that are usually employed for stability analysis. 

2 Ternary Mixture 

In this section we will first analyse the evaporation of a mixture of soil particles 
and water. As the liquid evaporates, both the temperature and the concen- 
tration change at the surface. This causes an instability to occur which drives 
convection. In order to analyse the occurrence of this stability, we first observe 
that the density of the liquid depends on both the temperature of this mixture 
and concentration of particles in water. However, in any real situation there 
will also be various solute impurities in the soil. Hence, we analyze a ternary 
mixture of water, soil particles and dissolved solutes. Thus, we can write, 

p = poll - a c (C - Co) + a' c (C - Q + a t (T - T )}, (1) 

where C is the concentration of soil particles and C is the concentration of 
solute impurities. The coefficients a c , a' c and at are the particle, solutal and 
thermal expansion coefficients, respectively, taken to be constants. We now let 
C —> -CAC + C , C -> -C'AC + C and T -> -TAT + T . We also define 
B c = a c AC,B' c = a' c AC and B t = a t AT and so we get, p = p [l + B C C + 
B' C C + B t T]. Our system is described by a divergenceless vector field, which 
represents the fluid velocity, d l Vi = 0. In order to write the momentum balance 
equation, it is useful to define a substantive derivative as follows 

Dvi = d t Vi + tP djVi, (2) 

where dt = d/dt, di = d/dx 1 , and d l di = d 2 . In continuum mechanics, this 
describes describes the time rate of change of some physical quantity for a 
material element subjected to a space-and-time-dependent velocity field. Now 
we can write the momentum balance equation as 

p Dvi = -d iP + d*ii[l - C/Cp}- 2 ^ + djVi) - pgXi. (3) 

Here p is the density, g is the acceleration due to gravity, p is a constant viscosity 
and C p is a maximum close packing concentration (shrinkage limit) |38j . The 
factor [1 — C/C p ]~ 2 denotes the dependence of viscosity on concentration. Along 
with this equation there are also the following equations, 

DT = nd 2 T, 

DC = d l d 11 [l-C/Cp]- 2 d l C + d 12 d 2 C', 
DC = d 21 2 C + d 22 d 2 C, (4) 

where k is the thermal diffusivity of the fluid, dn is the soil diffusion coef- 
ficient, d 22 is the solute diffusion coefficient, and d 2 \,d\ 2 are cross diffusion 
coefficients. We have neglected both the Dufour effect and Soret effect in the 
temperature equation as they are very weak relative to other effects. Here 



3 



the thermal diffusivity of the fluid, solute diffusion coefficient and the cross 
diffusion coefficients are constants. As we will be analysing the this system 
far away from C p , we will neglect the factor C/C p and so in this limit even 
soil diffusion coefficient and the viscosity is a constant. We will also use the 
Boussinesq approximation and set density constant in all terms except pg\. 
We can use the natural length scale I, which is the depth of the liquid, and 
thermal diffusivity k to non-dimensional these equations. To do so we will 
transform each quantity as x % — > lx\di — > l~ 1 di,Vi — > kl~ 1 Vi,t — > l 2 k~ x t, and 
p — > ppkl~ 2 + pog(l — AVi). The new velocity field again remains divergenceless, 
d l Vi = 0. The non-dimensional form of the momentum conservation equation 
and the continuity equation for temperature and concentration for constant 
viscosity and diffusion coefficients becomes 



Pr^Dvi = -d l p + d 2 v t + R c CX l +R' c C'\ l + R t TX t , 

DC = Le xl 1 d 2 C + Le^d 2 C, 

DC = Le^d 2 C + Le^d 2 C, 

DT = d 2 T, (5) 

where Pr = kpofi^ 1 is the Prandtl number, R c = p~ 1 k~ 1 B c pQgl 3 and R' c = 
/i _1 fc _1 B' c pogl 3 are the concentration Rayleigh numbers, Rt = p~ l k~ x B t pogl 3 
is the temperature Rayleigh number, Leu = fcd 1 ~ 1 1 ,£e22 = kd^ are the Lewis 
numbers and Le\2 = kd^ i -^ e 2i = kd%i are new Lewis number corresponding to 
cross diffusivity. We want to analyze the steady state version of these equations. 
In order to do so, we impose the following boundary conditions, v z = 0, d z v x = 
0, d z v v = 0, and C = C = T = 1 , (7(1) = C'(l) = T(l) = 0, p(l) = 0. We also 
take A; = (0, 0, 1) and = (x, y, z), thus our liquid is placed on a level plane. 
Under these boundary conditions we can obtain the steady state solutions. We 
will also add small perturbations to these steady state solutions. Now the steady 
state solutions with perturbations added to them is are given by 

Vi = u, 

c = i-\ i r i + e c 



C = 1 - A 
T = 1 - A 



(R c + R' c + R t ) , _ .j |2 
2 

Now we define D^ = a~ 1 dt — b~ 1 d 2 , and so we get 



V = - ^^^ |1-AV,1 2 + ^. (6) 



D( r Ui + di<t>- RcOXi = 0, 

D\ eXX B c - Le^d 2 9' c - = 0, 

D 1 Le22 e' c -Le^d 2 e c -\ l u l = 0, 

D\e t -Xui = 0. (7) 

We define an operator E which takes the curl of any vector field a % twice, 

Ea? = e mn e nkl d m d k a l 

= (5iSr-6j6Z)d m d k a l 

= d l d p a p ~d 2 a l . (8) 
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So, the gradient of a scalar field vanishes when its curl is taken twice, Ed l <p = 
0, because in the expression Ed l <f> = e lmn e n kid m d k d l 4> there is a contraction 
between a pair of symmetric and antisymmetric tensor indices. We also have 
Eui = —d 2 u.i, because Ui is divergenceless, and X l (EXi(j)) = \ l di\ % dj<j) — d 2 4> = 
d 2 <f>. Now acting on the momentum balance by E and contracting it with Ai, 
we get 

Df r d 2 \ l Ul - R c d 2 9 c - R' c d 2 9' c - R t d 2 9 t = 0. (9) 
Now we define the following differential operator 

V = Dl e22 Dl ell -Le^Le^d\ 
Vi = (Dl^-Le^d 2 ), 



2>2 = ODLii-W^ 2 ), 



Now we have 



So, we have 



because 



(10) 



and the following Green's function's 

V G c (r,r') = 6 3 (r-r'), 

D\G t {r,r') = 5 3 {r~r'). (11) 



V 9 c = V 1 Xu. l , 
V 8' c = V 2 X l u u 

D\0 t = A*«i. (12) 



c (r) = J d 3 r'G c (r,r , )V 1 X i u i (7 J ), 
0' c (r) = J d 3 r'G c (r,r')2? 2 A 4 Ml (r'), 
6 t (r) = J d 3 r / Gf t (r,r / )A i Ui(r). (13) 

T> e c (r) = f dr / 2> G e (r,r')2) 1 A < u i (r'), 
d 3 r'6 3 (r - r')£>iAX(r') 

V 6 c (r) = f dr'P G c (r,r')^ 2 AV t (r'), 

d 3 r'S 3 {r - r')V 2 \ l u t {r') 
= V 2 X l u t (r), 
D\9 t {r) = [ dr , D{Gt(ry)\ i m(r f ) 

d 3 r'5 3 (r -r')A 4 w 4 (r') 
\ % Ui(r). (14) 
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Thus, we can write, 

D^ r d 2 X i u i (r) = R c d 2 j d 3 r'G c (r,r')I'iA i u i (r') 
+R' c d 2 J dVG c (r,r')P 2 AV 2 (r') 
+R t d 2 J (Pr'Gtiry^Uiir'). (15) 

Multiplying by Vq and D{, we get 

T> DlD[ r d 2 X i u i = i? c a 2 J DiPiA 4 u 4 + ^9 2 L> 1 1 2? 2 A t M J 

+i? t 9 2 P A l Ui . (16) 

Now using the boundary conditions that the even derivatives of X l di vanishes 
on X l Ui, we write the solution as [57] . 

Xiii = A sin rnr exp[i(k x x + k y ) + at], (17) 

becomes and we also define 

fc 2 = fc 2 + fc 2 , k 2 n = k 2 + n 2 TT 2 . (18) 

Thus the characteristic equation is 

C{k m ,k,a)=0, (19) 

where 

C(k m ,k, cr) = ((a + k^Leg^a + fc 2 ^^ 1 ) - LegLegk^) 
xia + k^aPr- 1 + k 2 n )k 2 n 
-R c (a + k 2 n ){{a + k 2 n Leg) - Legk 2 n )k 2 
-R' c (a + k 2 n )((a + k 2 n Leg) - Legk 2 n )k 2 
-R t ((a + k 2 n Leg){<7 + k 2 n Leg) 
-LegLe^k^k 2 . (20) 

Assuming the principle of exchange of stabilities critical behavior is obtained by 
setting cr = [37], we have 

C{k m ,k,0) = 0, (21) 

where 

C(k m , k, 0) = {Leg Leg - Leg Leg)k 6 n k~ 2 
—R c (Le 22 — Le 12 ) 
-R' c (Leg - Leg) 

-R t {LegLeg - LegLeg), (22) 

Now we can define a effective Raylcigh number as 

R = [R c (Leg - Leg) + R' c (Leg - Leg) 
+R t (LegLeg - LegLeg)] 

x(LegLeg - LegLeg)- 1 . (23) 
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Hence, we can write R = k^k 2 . The lowest value is for n = 1 and the instability 
starts at 

This gives us A: 2 = 7r 2 /2, and the corresponding value of R will be given by 
R = 657.5 - 10 3 . 

This is when the instability will start. This convection in ternary mixtures 
can be used to analyse interesting geological phenomenon. This is because soils 
can exhibit semi-permeability and osmosis through these semipermeable soils 
can be driven by gradient in salt concentration [39j 0D1 [41] . In fact, an excess 
pressure has been observed to exist during the diffusion of salty water through 
certain soils [121113]. Thus, it will be interesting to analyse this phenomenon us- 
ing stability analysis for ternary systems. In the limit when there is no cross dif- 
fusion, Le^2 = Le-2i = 0, we have the expected result R — R c L\i + R' C L22 + Rt- 
Thus, in this limit, the temperature Raylcigh number and all the concentration 
Rayleigh multiplied by there respective Lewis numbers add up to give the ef- 
fective Rayleigh number. So, the effective Rayleigh number is enhanced by the 
existence of a salt gradient. This, also implies that convection can start much 
earlier in salty water than pure water. Thus, the evaporation of salty water at 
the surface of soils can set up interesting instabilities, governed by the ternary 
characteristic equation. 



3 Binary Mixture 

In this section we will neglect the effect due to impurities. However, we are not 
able to directly set Le^ = ^ e 2i = ^ e i2 = 0' m the effective Rayleigh num- 
ber. This is because the Green's function used to derive this effective Raylcigh 
number is not well defined for these values. This is because the differential op- 
erator T>o has a zero eigenvalue for these values of the Lewis numbers. So, its 
inverse does not exist. Thus, we can not set Le^ = Le^i = Le^ 2 l = 0, in the 
effective Rayleigh number obtained by using this Green's function. To obtain a 
correct result we will have to repeat the above analysis for the binary solutions. 
However, it may be noted that all the analysis of the previous section, except 
the derivation of the Green's function remains, remains well defined, if we set 
^ e 22 = = Le^ = 0. So, now we set, Leu = Le, Le^2 = L e 2i = ^ e Y2 = 

0, p = Po[l + B C C + B t T], and use the following equations 

Pr- l Dv % = -d lP + d 2 v l +R c C\ l + R t TX l , 
DC = Le~ 1 d 2 C, 

DT = d 2 T. (25) 

If we repeat the above analysis, we will obtain the following pcrturbative equa- 
tions 

D Pr d 2 \ l Ul - R c d 2 9 c - R t d 2 6 t = 0, 
D\ e 9 c -\ l u t = 0, 
D\6 t - Km = 0. 

(26) 
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Now we again define the following Green's functions 

D l Le G(rS) = 5 3 (r-r>), 

D\G t {r,r') = S 3 (r-r'). (27) 

Thus, we can write 

c {r) = J d 3 r'G{r,r')y Ul (r'), 

t (r) = J dVGf t (r, OAVrO, ( 28 ) 

because 



DlMr) = I <Pr'D 1 Le G(r,r')\ i u i (r') 
d 3 r'5{r -r')\ l Ul {r') 

D\e t {r) = I ' d 3 r , D{G t (r,r')\ i u i {r') 

d 3 r'5{r -r')\ l Ul {r') 
X l u t (r). (29) 



We can write 



D Pr d 2 X i u i (r) = R c d 2 J d 3 r'G{r,r')X l Ul (r') 

+R t 3 2 f d 3 r'G t (r,r')X l u t (r'). (30) 



Now acting on this equation by D\ and D\ e , we get 

DlD 1 Le D[ r d 2 X i u i = R c d 2 D\\ l u % + R t d 2 Dl e \ l u t {r'). (31) 

So, again using the boundary conditions that the even derivatives of \ l di van- 
ishes on X l Ui, we write the solution as [57] 

\ l Ui = A sin rnr exp[i(k x x + k y y) + at], (32) 

Thus, we get 

C(k„,k,a) =0, (33) 

where 

C(k n ,k,a) = (a + klLe-^iPr^a + kDia + kX 

-(a + k 2 n Le- v )k 2 R t - (a + kl)k 2 R c . (34) 

So, in the case of temperature dependence critical behavior is obtained by setting 
(7 = 0, we have [57] 

C(fc m ,fc,0) = 0, (35) 



where 



C(kn,k,0) = k®k~ 2 Le~ 1 — Le~ l R t - R c - (36) 

Now here the effective Rayleigh number is 

R= R t + LeR c (37) 

Hence, we can write R — k%kr 2 . The lowest value is for n = 1 and the instability 
starts at 

r) /? 

aP = - ™ 
This gives us k 2 = 7r 2 /2, and the corresponding value of R will be again be 
given by R = 657.5 ~ 10 3 . 

Now the actual value for a binary mixture of bentonite and water are of the 
order, p ~ 10 3 , fj, ~ 10~ 3 , k ~ 1CT 7 . Also the depth of the liquid I ~ 1CT 2 [22] 
and B c ~ l,B t ~ 10 _1 [H], so we have R t ~ 10 3 , and R c ~ 10 6 . However, in 
the stability analysis the product of R c and ie contributes to the occurrence 
of the instability and Le ~ 10 2 . So, the contribution from the concentration 
gradient R c Le ~ 10 s is much more than temperature gradient. The Marangoni 
effect is also measured by the Marangoni number M ~ 10 3 [45] . Hence, the 
concentration gradient is the most dominant factor for convection to occur. 
The calculated value of the effective Rayleigh number due to the concentration 
gradient is much larger than the theoretical limit for the convection to occur. 
Hence, we hypothesize that it is responsible for the instability in the drying of 
a binary mixture of bentonite and water. 

It has been observed that in a binary mixture of bentonite and water initially 
the convection takes place on the surface, but, eventually a layer forms on the 
surface and this inhibited further convection from taking place on the surface 
[22] . However, the convection continues in the lower depths of the soil. This ob- 
servation rules out the possibility that this convection is due to the Marangoni 
effect. This is because in that case it would be surface driven phenomenon and 
would not continue at depth after halting at the surface. It will be interesting 
to explore the possibility that this behavior is due to the fact that in a soil the 
diffusivity and viscosity are strong functions of particle concentration. That is, 
the effective viscosity and diffusivity of a soil slurry increase dramatically with 
increasing particle concentration, and thus at sufficiently high particle concen- 
trations the effective Rayleigh number will be reduced to below the critical value 
necessary for convection. 



4 Conclusion 

In this paper we have analysed the drying of a ternary mixture of soil particles 
and water mixed with impurities. We have demonstrated that the concentra- 
tion gradient that forms during the process of evaporation of this mixture is 
the main driving force for the occurrence of a Rayleigh-Benard instability. This 
instability causes convection to take place, which in turn may have significant 
implications for the drying of soil crusts and the formation of patterned ground. 
In the ternary case the use of Green's functions allowed for a straight-forward 
calculation of the characteristic equation. For the binary case a separate calcu- 
lation was required. This is because the differential operator that was inverted 
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using the Green's function in the ternary case, yields a zero eigenvalue in the 
binary case. Thus, it cannot be inverted and the ternary Green's function be- 
comes ill defined for the binary case. As noted in the introduction, soil loss in 
agricultural areas depends critically on the formation of surface seals and crusts 
and the subsequent cracking of these crusts during drying of the ground. It is 
possible that the convection patterns give rise to weak points, which on drying 
gives rise to cracks [22] . Furthermore, patterns observed during the freezing 
and thawing of ground around the Arctic circle have also been explained as 
owing to convection of the soil. In this paper it was proposed that evaporation 
from the surface of thawing ground causes an unstable density profile to develop 
potentially leading to convection. The occurrence of similar patterns on Mars 
suggests that in the past the temperature on Mars may have been large enough 
for the ice to melt and evaporation to take place. It may be noted that it has 
been observed that drying on slanting slopes gives rise to rolls and drying on flat 
surface gives rise to hexagonal structures |32j . It will be interesting to perform 
stability analysis for these specific physical situations and see if the formation 
of these particular patterns can be explained. 
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